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, We propose a multipole representation of the Fermi-Dirac function and the Fermi operator, and 

use this representation to develop algorithms for electronic structure analysis of metallic systems. 
' The new algorithm is quite simple and efficient. Its computational cost scales logarithmically with 

, /3Ae where /3 is the inverse temperature, and Ae is the width of the spectrum of the discretized 

. Hamiltonian matrix. 
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^ ; I. INTRODUCTION 

The Fermi operator, i.e. the Fermi-Dirac fmiction of the system Hamiltonian, is a fundamental quantity in the 
quantum mechanics of many-electron systems and is ubiquitous in condensed matter physics. In the last decade the 
development of accurate and numerically efficient representations of the Fermi operator has attracted lot of attention 
. in the quest for linear scaling electronic structure methods based on effective one-electron Hamiltonians. These 
approaches have numerical cost that scales linearly with N, the number of electrons, and thus hold the promise of 
making quantum mechanical calculations of large systems feasible. Achieving linear scaling in realistic calculations 
■ O ' is very challenging. Formulations based on the Fermi operator are appealing because this operator gives directly the 
C ] single particle density matrix without the need for Hamiltonian diagonalization. At finite temperature the density 
matrix can be expanded in terms of finite powers of the Hamiltonian, requiring computations that scale linearly with 
A'^ owing to the sparse character of the Hamiltonian matrix [Tsj. These properties of the Fermi operator are valid not 
only for insulators but also for metals, making formulations based on the Fermi operator particularly attractive. 
' Electronic structure algorithms usinga Fermi operator expansion (FOE) were introduced by Baroni and Giannozzi 
, [H and Goedecker and co-workers 0, 13| (see also the review article Q). These authors proposed polynomial and 
' rational approximations of the Fermi operator. Major improvements were made recently in a series of publications by 
, Parrinello and co-authors 0, 0, 0, H H fioj , in which a new form of Fermi operator expansion was introduced based 
• on the grand canonical formalism. 

[ From the viewpoint of efficiency, a major concern is the cost of the representations of the Fermi operator as a 
I • function of /3Ae where /? is the inverse temperature and Ae is the spectral width of the Hamiltonian matrix. The cost 
QQ ■ of the original FOE proposed by Goedecker et al scales as /?Ae. The fast polynomial summation technique introduced 
, by Head-Gordon et aZ[ll|, reduces the cost to (/3Ae)^/^. The cost of the hybrid algorithm proposed by Parrinello 
' et al in a, recent preprint scales as (/3Ae)^/'^. 
. ^ , The main purpose of this paper is to present a strategy that reduces the cost to logarithmic scaling (In (3Ae) ■ 
■ (lnln/3Ae), thus greatly improving the efficiency and accuracy of numerical FOEs. Our approach is based on the 
J-j [ exact pole expansion of the Fermi-Dirac function which underlies the Matsubara formalism of finite temperature 
. . . 1 Green's functions in many-body physics [l3| . It is natural to consider a multipole representation of this expansion to 
achieve better efficiency, as was done in the fast multipole method Indeed, as we will show below, the multipole 
expansion that we propose does achieve logarithmic scaling. We believe that this representation will be quite useful 
both as a theoretical tool and as a starting point for computations. As an application of the new formalism, we 
present an algorithm for electronic structure calculation that has the potential to become an efficient linear scaling 
algorithm for metallic systems. 

The remaining of the paper is organized as follows. In the next section, we introduce the multipole representation 
of Fermi operator. In Section lllli we present the FOE algorithm based on the multipole representation and analyze 
its computational cost. Three examples illustrating the algorithm are discussed in Section [TVl We conclude the paper 
with some remarks on future directions. 
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II. MULTIPOLE REPRESENTATION FOR THE FERMI OPERATOR 



Given the effective one-particle Hamiltonian H, the inverse temperature /3 = and the chemical potential ji, 

the finite temperature single-particle density matrix of the system is given by the Fermi operator 



2 f P 

7— 7r-= = 1 — tanhl — (H — u) ) , 

l4-exp(/3(H-/^)) V2^ 



(1) 



where tanh is the hypcrbohc tangent function. The Matsubara representation of the Fermi-Dirac function is given by 

00 ^ 

''-'-'"^-i: ^(H-rt-(2--i).. - 

This representation originates from the pole expansion (see, for example, (isl. [l6t) of the meromorphic function tanh 

^ 1 



tanh z = 



2/-1, 



;=-oo " 2 

In particular, for x real (which is the case when H is self-adjoint), we have 



00 _^ 

tanh X = 4Re >^ 

^ 2x - (2; - 1 

1=1 ^ 



(3) 



(4) 



To make the paper self-contained, we provide in Appendix A a simple derivation of this representation. It should be 
emphasized that ([2]) is exact. We notice that the expansion ([2]) can also be understood as the P — > 00 limit of an 
exact Fermi operator expansion proposed by Parrinello and co-authors in 0, 0, 0, H, 0, ■ 

The summation in Q can be seen as a summation of residues contributed from the poles {{21 — l)7ri}, with I a 
positive integer, on the imaginary axis. This suggests to look for a multipole expansion of the contributions from 
the poles, as done in the fast multipole method (FMM) [l^. To do so, we use a dyadic decomposition of the poles, 
in which the n-th group contains terms from I = 2"~^ to Z = 2" — 1, for a total of 2"~^ terms (see Figure [T] for 
illustration). We decompose the summation in Eq.Q accordingly, denoting x = /3(H — fi) for simplicity 



y 

'p^x- {21- l)TTi 



00 2" 

E E 

n=l Z=2' 



X- {21- l)7:i 



71=1 



(5) 



The basic idea is to combine the simple poles into a set of multipoles at / = where In is taken as the midpoint of 
the interval [2"-i,2" - 1] 



3 • 2"-i - 1 



(6) 



Then the S'„ term in the above equation can be written as 

2"-l 



E 

Z=2"- 
2"-l 

E 

=2"- 

2"-: 

E 



1 



1=2 

2"-l 



X — {2ln — l)7ri — 2{l — ln)TTi 

1 . 2{l- 

Diri ^\ 



^ X - {2ln - l)Tri ^^^x - {2ln - l)m' 

1 2{l-lr,)TTi Y 

„-l)7^^ ^U-(2/„-l)7rJ 



(7) 



X - {2ln - 1)..- 
2"-l 



E : 



1 



;=2' 



2{l — Z„)7ri 

1^ X- {21- \)m \x- {2ln - 1)7 



In deriving Eq. ([7]) we used the result for the sum of a geometric series. Using the fact that x is real, the second term 
in Eq. ([7]) can be bounded by 



E 



^_^ \ X — {21 — X — {2ln — l)7rz 



2{l — ^„)7ri 



-^E_J(2/-l)7r| 



2{l-ln) P ^ 

2L-1 - 2ti iP 



(8) 
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Therefore, we can approximate the sum Sn by the first P terms, and the error decays exponentially with P 

2"-l P-1 



2{l — ln)T^i 



^»(^) E X - (2L - Dttz E 



uniformly in x. The above analysis is of course standard from the view point of the fast multipole method [14||. The 
overall philosophy is also similar: given a preset error tolerance, one selects P, the number of terms to retain in Sn, 
according to Eq. 

Interestingly, the remainder of the sum in Eq. ^ from I — m to oo has an explicit expression 

Re ^ — = — Im V I "I - - + -x] , (10) 

^ 2x- (21 - l)iTT 27r ^ V 2 TT / ' ^ ' 

l—m 

where ip is the digamma function V'(-z) = r'(z)/r(z). It is well known [l^ that the digamma function has the following 
asymptotic expansion 

^{z)^\n{z)-^^-^ + 0(\), |argz| <^ and |z| ^oo. (11) 



Therefore, 



2z \2z'^ 



1 4 \ . . / 1 i \ „/ 1 



Im th [ m \ X ^ Im In rii x ] + O , 

^' 2 n J \ 2 TT J Vj7l2 

"'^^^""( (2J-l)J +^(;^ 



(12) 



Figure [5] shows that the asymptotic approximation (|12p is already rather accurate when m = 10. 

Eq. p2p also shows the effectiveness of the multipole representation from the viewpoint of traditional polynomial 
approximations. At zero temperature, the Fermi-Dirac function is a step function that cannot be accurately approx- 
imated by any finite order polynomial. At finite but low temperature, it is a continuous function with a very large 
derivative at x = 0, i.e. when the energy equals the chemical potential ^. The magnitude of this derivative becomes 
smaller and, correspondingly, the Fermi function becomes smoother as the temperature is raised. One can use the 
value of the derivative of the Fermi function at a; = to measure the difficulty of an FOE. After eliminating the first m 
terms in the expansion, Eq. p2p shows that asymptotically the derivative is multiplied by the factor ^2rn-i)-n ' ^^ich 

is equivalent to a rescaling of the temperature by the same factor. In particular, if we explicitly include the first 2^ 
terms in the multipole representation of the Fermi operator, we are left with a remainder which is well approximated 
by Eq. (fT2l) . so that, effectively, the difficulty is reduced by a factor 2^ . As a matter of fact standard polynomials 
approximations, such as the Chebyshev expansion, can be used to efficiently represent the remainder in Eq. (jlOp even 
at very low temperature. 

In summary, we arrive at the following multipole representation for the Fermi operator 

N 2'^-! ^ P-1 nfl 1 \ ■ 

1 ■r-^ / 2[l — ln)Tn 



p=l-4Re^^ 



^Im (2^ - i + ^/3(H - /i)) + 0(iV/3^). (13) 

The multipole part is evaluated directly as discussed below, and the remainder is evaluated with the standard poly- 
nomial method. 



III. NUMERICAL CALCULATION AND ERROR ANALYSIS 



To show the power of the multipole expansion, we discuss a possible algorithm to compute the Fermi operator in 
electronic structure calculations and present a detailed analysis of its cost in terms of /3Ae. Given the Hamiltonian 
matrix H, it is straightforward to compute the density matrix p from the multipole expansion if we can calculate the 
Green's functions = [/3(H — //) — (2Z„ — \)t:i\~^ for different n. 
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A possible way to calculate the inverse matrices is by the Newton-Schulz iteration. For any non-degenerate matrix 
A, the Newton-Schulz iteration computes the inverse B = A~^ as 

Bfc+i = 2Bfe - BfcABfc. (14) 

The iteration error is measured by the spectral radius, i.e. the eigenvalue of largest magnitude, of the matrix I — BBfc 
where I is the identity matrix. In the following we denote the spectral radius of the matrix A by cr(A). Then the 
spectral radius at the k-th step of the Newton-Schulz iteration is Rfe = a{l — BBfc) and 

a(Rk+i) = a(Rfe)' = a{Rof^\ (15) 

Thus the Newton-Schulz iteration has quadratic convergence. With a proper choice of the initial guess (see 0), the 
number of iterations required to converge is bounded by a constant, and this constant depends only on the target 
accuracy. 

The remainder, i.e. the term associated to the digamma function in Eq. (|13p. can be evaluated by standard 
polynomial approximations such as the Chebyshev expansion. The order of Chebyshev polynomials needed for a 
given target accuracy is proportional to /3Ae/2^+^ (see [13, Appendix]). 

Except for the error coming from the truncated multipole representation, the main source of error in applications 
comes from the numerical approximation of the Green's functions B/^ . To understand the impact of this numerical 
error on the representation of the Fermi operator, let us rewrite 

2"-l P-l P-l 2"-l 

S„= B/„^(-2(Z-;„)7rzB,J'' = ^Br+i ^ (-2(Z - ;„)7rz)^ 

The factor '^i{—'2{l — ln)'r:iY is large, but we can control the total error in Sn in terms of the spectral radius 
(t(B/,^ — B;^). Here B;^ is the numerical estimate of B;^. 
The error is bounded by 

p-l p-l 
ct(S„ - S„) < ^ 2"-1(2"-V)''ct(B''+i - b"^') < ^(2"-V)''+V(B"+i - b"^'), (16) 

v=0 i/=0 

where we have omitted the subscript In in B/^ and in B;^. In what follows the quantity 'Yli^Zo{'^^^~^''^Y^^'^{^^^^ ~ 
^i^+i, 

B j will be denoted by ep. Then we have 

p-l 



ep = ^(2"-V)''+V((B'' - b'')B + (b'' - B'^)(B - B) + B''(B - B)) 

i/=0 

p-l p-l 
< (2"- V)''+i ((t(B) + (t(B - B))cr(B^ - b") ^ (2"-i7r)^+V(B)''cr(B - B). 



u=l i/=0 

Here we took into account the fact that the v ~ Q term in the first summation is equal to zero, and have used the 
properties (t(A H- B) < cr(A) + cr(B), and o-(AB) < cr(A)cr(B), respectively. 

Noting that 2"^^7ro'(B;^ ) < 1/3 and changing v io v + 1 in the first summation, we can rewrite ep as 

p-2 p-l 
ep < (r + 2"-i7ra(B - B)) ^ (2"-i7r)'^+i(T(B'^+i - b''^^) + ^ —T-^Tia{^ - B) 

< (i + 2"-i7ra(B - B))ep_i + -2"-i7rCT(B - B) ^^^^ 
3 2 

,1 , 3 
= (g + ei)ep_i + -ei. 

In the last equality, we used the fact that ei ~ 2"^^7ra(B — B). Therefore, the error ep satisfies the following recursion 
formula 

3ei/2 ,1 , / 3ei/2 \ A ,p . f 3ei/2 \ 
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Taking ei < |, we have 



ep < ei = 2"-Vcr(B -B). (20) 



Therefore, using Eq. (|15p wc find that the number k of Newton-Schulz iterations must be bounded as dictated by 
the following inequality in order for the error (t(S„ — S„) to be < 10~^ /N. 

2"-V(Ro)2^ < (21) 

Here we have used the fact that cr(B;^) < I/tt for any n. Each Newton-Schulz iteration requires two matrix by matrix 
multiplications, and the number of matrix by matrix multiplications needed in the Newton-Schulz iteration for B;^ 
with n < is bounded by 



log2 cr (Ro) 

To obtain a target accuracy (j{p — p) < 10^^ for a numerical estimate p of the density matrix, taking into account 
the operational cost of calculating the remainder and the direct multipole summation in the FOE, the number of 
matrix by matrix multiplications hmm is bounded by 

riMM < 2^1og2 N + CiN + C22-^-i/3Ae. (23) 

Here we used the property: log2(a; + y) < log2 x + log2 y when x,y > 2, and defined the constant Ci as follows: 



n=l \ 



log2Cr((Ro)/„' 



The dependence on 2~^~-^PAe in the last term on the right hand side of ((23|) comes from Chebyshev expansion used 
to calculate the remainder. From numerical calculations on model systems, the constant Ci and C2 will be shown to 
be rather small. Finally, choosing N cx ln(/3Ae), we obtain 

nuM oc (In /3Ae) • (In In (3Ae) (25) 

with a small prefactor. 



IV. NUMERICAL EXAMPLES 



We illustrate the algorithm in three simple cases. The first is an off-lattice one dimensional model defined in a 
supercell with periodic boundary conditions. In this example, we discretize the Hamiltonian with the finite difference 
method, resulting in a very broad spectrum with a width of about 2000eV, and we choose a temperature as low as 
32K. In the second example wc consider a nearest neighbor tight binding Hamiltonian in a three dimensional simple 
cubic lattice and set the temperature to lOOK. In the third example we consider a three dimensional Anderson model 
with random on-site energy on a simple cubic lattice at lOOK. 



A. One dimensional model with large spectral width 



In this example, a one dimensional crystal is described by a periodic supercell with 10 atoms, evenly spaced. We 
take the distance between adjacent atoms to be a = 5.29A. The one-particle Hamiltonian is given by 



1 92 



The potential V is given by a sum of Gaussians centered at the atoms with width a = 1.32A and depth Vo = 13.6eV. 
The kinetic energy is discretized using a simple 3-point finite difference formula, resulting in a Hamiltonian H with a 
discrete eigenvalue spectrum with lower and upper eigenvalues equal to e_ — 6.76eV and e+ = 1959eV, respectively. 
Various temperatures from 1024K to 32K were tried. Figure [3] reports the linear-log graph of tt-mm, the number of 
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matrix by matrix multiplications needed to evaluate the density matrix using our FOE. versus /?Ae, with /3Ae plotted 
in a logarithmic scale. The logarithmic dependence can be clearly seen. The prefactor of the logarithmic dependence 
is rather small: when /3Ae is doubled, a number of additional matrix multiplications equal to 17 is required to achieve 
two-digit accuracy {D = 2), a number equal to 19 is needed for £> = 4, and a number equal to 21 is needed for D = 6, 
respectively. The observed L'-dependence of the number of matrix multiplications agrees well with the prediction in 
(1231). 

In order to assess the validity of the criterion for the number of matrix multiplications given in Eq. (23), we report 
in Table [J the calculated relative energy error and relative density error, respectively, at different temperatures, when 
the number of matrix multiplications is bounded as in formula (23) using different values for D. The relative energy 
error, Aeid, measures the accuracy in the calculation of the total electronic energy corresponding to the supercell 
E = Tr(pH). It is defined as 

A.„, = E^. ,27, 

Similarly the relative error in the density function in real space is defined as 

Ap.e, = (28) 
Trp 

Because Trp = N^, where Ne is the total number of electrons in the supercell, Aprei is the same as the density 
error per electron. Table |T] shows that for all the values of f3Ae, our algorithm gives a numerical accuracy that is even 
better than the target accuracy D. This is not surprising because our theoretical analysis was based on the most 
conservative error estimates. 



B. Periodic three dimensional tight-binding model 



In this example we consider a periodic three dimensional single-band tight-binding Hamiltonian in a simple cubic 
lattice. The Hamiltonian, which can be viewed as the discrctized form of a free-particle Hamiltonian, is given in 
second quantized notation by: 

H = -t^c+c„ (29) 

<hj> 

where the sum includes the nearest neighbors only. Choosing a value of 2.27eV for the hopping parameter t the band 
extrema occur at e_|_ = 13.606cV, and at e_ ~ — 13.606eV, respectively. In the numerical calculation we consider a 
periodically repeated supercell with 1000 sites and chose a value of lOOK for the temperature. Table [Til shows the 
dependence of umm, Aerci, and Aprci on the chemical potential ^, for different D choices. Compared to the previous 
one dimensional example in which (3Ae was as large as 7.12 x 10^, here /?Ae = 1600 due to the much smaller spectral 
width of the tight-binding Hamiltonian. When fi = the chemical potential lies exactly in the middle of the spectrum. 
This symmetry leads to a relative error as low as 10~^^ for the density function. 



C. Three dimensional disordered Anderson model 



In this example we consider an Anderson model with on-site disorder on a simple cubic lattice. The Hamiltonian 
is given by 

H=-i ^ c,+c, +5]e,c+c,. (30) 

<ij> i 

This Hamiltonian contains random on-site energies ei uniformly distributed in the interval [— 1.13eV, 1.13eV], and we 
use the same hopping parameter t as in the previous (ordered) example. In the numerical calculation we consider, as 
before, a supercell with 1000 sites with periodic boundary conditions, and choose again a temperature of lOOK. In 
one realization of disorder corresponding to a particular set of random on-site energies, the spectrum has extrema at 
e+ = 13.619eV and at e_ = — 13.676eV. The effect of disorder on the density function is remarkable: while in the 
periodic tight-binding case the density was uniform, having the same constant value at all the lattice sites, now the 
density is a random function in the lattice sites within the supercell. Table IIIII reports for the disordered model the 
same data that were reported in Table |ll] for the ordered model. We see that the accuracy of our numerical FOE is 
the same in the two cases, irrespective of disorder. The only difference is that the super convergence due to symmetry 
for ^ = no longer exists in the disordered case. 
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V. CONCLUSION 



We proposed a multipole representation for the Fermi operator. Based on this expansion, a rather simple and 
efficient algorithm was developed for electronic structure analysis. We have shown that the number of number of 
matrix by matrix multiplication that are needed scales as (ln/?Ae) • (lnln/3Ae) with very small overhead. Numerical 
examples show that the algorithm is promising and has the potential to be applied to metallic systems. 

We have only considered the number of matrix matrix multiplications as a measure for the computational cost. The 
real operational count should of course take into account the cost of multiplying two matrices, and hence depends on 
how the matrices are represented. This is work in progress. 



APPENDIX A: MITTAG-LEFFLER'S THEOREM AND POLE EXPANSION FOR HYPERBOLIC 

TANGENT FUNCTION 

To obtain the pole expansion for hyperbolic tangent function tanh(z) , we need a sp ecial case of the general Mittag- 
Leffler's theorem on the expansions of meromorphic functions (see, for example, |l5l. [l6j). 

Theorem 1 (Mittag-LefFler). // a function f{z) analytic at the origin has no singularities other than poles for finite 
z, and if we can choose a sequence of contours Cm about z = tending to infinity, such that \ f{z)\ < M on Cm and 
Idz/^l is uniformly bounded, then we have 

f{z) = /(O) + lim {Pmiz) - PmiO)}, (Al) 

m — ''Oo 

where Pm{z) is the sum of the principal parts of f{z) at all poles within Cm- 

For tanh(2) = -, it is analytic at the origin and tanh(O) = 0. The function has simple poles at 

exp(z) + exp(— z) 

zi = [l — l/2)7ri, I G Z with principle parts {z ~ zi)~^. Let us take the contours as 

Cm = {x ± imn \ \x\ < mn} U {±mn + iy \ \y\ < mn}, m G Z+, 
it is then easy to verify that Cm satisfy the conditions in Theorem [1] According to Theorem [1] 

™ 1 1 

tanh(z) = tanh(O) + lim V f h — V (A2) 

l=-m+l ' ' 

By symmetry of zi, the second term within the brackets cancels, and we arrive at 
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FIG. 1: (color online) Illustration of the pole decomposition (|13|l . From 2" to 2"^^ — 1 poles are grouped together as shown 
the figure. The spectrum is indicated by the red line on the real axis. 
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FIG. 2: (color online) The function Im tp (m — i + ■^a::)(red circle), i.e. the remainder of the pole expansion in Eq. (|13|) 
compared with the function arctan ^ (2rn-i)n ) (blue solid line) for m — W 
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FIG. 3: (color online) Linear-log plot of the number of matrix matrix multiplications umm versus /3Ae. njviM depends logarith- 
mically on pAe with a small constant prefactor. 
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D = 4 
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W24:K 

512K 

256K 

128K 

6AK 

32K 


2.22 X lO* 
4.44 X 10* 
8.89 X 10* 
1.78 X 10^ 
3.56 X lO'^ 
7.12 X 10^ 


1.64 X lO"'' 

1.73 X 10"^ 
1.78 X 10"^ 

1.74 X 10"^ 

1.75 X 10"^ 

1.76 X 10"^ 


5.98 X 10"'= 
6.49 X 10"'^ 
6.83 X IQ-*^ 
6.55 X IQ-^ 
6.62 X IQ-'^ 
6.66 X 10"^ 


3.31 X IQ-'' 
3.70 X IQ-* 
3.96 X IQ-* 
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2.52 X 10"^ 
2.62 X 10"'' 
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2.70 X 10"** 
2.64 X lO"'' 


1.50 X lO"*" 
1.74 X 10"® 
1.81 X 10"® 
1.95 X 10"® 
1.86 X 10"® 
1.80 X 10"® 



TABLE I: One dimensional Hamiltonian model of Section FlV Al Relative energy error Aeioi and relative density error Aproi 
for a large range of values of /3Ae and several values of D. 
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-13 


2.52 X 10" 


15 


O.OOeV 


305 


1.55 X 10 


-9 


6.26 X 10" 


19 


357 


2.98 X 10" 


-15 


6.26 X 10" 


19 


5.44eV 


308 


1.45 X 10 


-8 


1.34 X 10" 


12 


356 


5.36 X 10" 


-13 


1.07 X 10" 


16 


10.88eV 


320 


1.69 X 10 


-8 


1.78 X 10" 


13 


376 


1.09 X 10" 


-12 


1.80 X 10" 


17 



TABLE II: Three dimensional periodic tight binding model of Section [IV Bl Number of matrix matrix multiplications tt-mm, 
relative energy error Aeroi, and relative density error Aproi. For p = 0, the algorithm achieves machine accuracy for the 
absolute error of the density function as a consequence of symmetry. 







D 


= 4 








D 


= 8 




riMM 


Aercl 


Aprcl 


JlMM 


Aercl 


Aprcl 


-10.88eV 


320 


5.16 X 10 


-9 


1.72 X 10" 


lU 


376 


3.16 X 10" 


-13 


2.59 X 10" 


14 


-5.44eV 


308 


4.75 X 10 


-9 


2.43 X 10" 


11 


356 


3.71 X 10" 


-13 


1.48 X 10" 


15 


O.OOeV 


305 


8.08 X 10" 


10 


9.50 X 10" 


13 


357 


1.76 X 10" 


-14 


2.39 X 10" 


17 


5.44eV 


308 


1.01 X 10 


-8 


1.22 X 10" 


12 


356 


3.57 X 10 


-13 


8.05 X 10" 


17 


10.88eV 


320 


1.30 X 10 


-8 


1.56 X 10" 


13 


376 


9.56 X 10" 


-13 


1.83 X 10" 


17 



TABLE III: Three dimensional Anderson model with on-site disorder discussed in Section II V CI Number of matrix matrix 
multiplications timm, relative energy error Aerd, and relative density error Aprd. 



